#ifndef Algebra_ODESystem_h
#define Algebra_ODESystem_h

#include "AlgebraModule.h"
#include "AlgebraicTypes.h"

SG_NAMESPACE_OPEN

namespace Algebra
{
    /**
     * \addtogroup ODE
     * @{
     */

    /**
     * @brief Represent the first-order ODEs \f$ \dot{u} = G\left ( t, u \right ) \f$.
     *
     * To go further, DAEs can be viewed as a generalization of an ordinary differential equations (ODEs).
     * A differential-algebraic equation (DAE) is an equation involving an unknown function and
     * its derivatives.
     *
     * \f[
     *  \mathbf{f} (t, u,\frac{du}{dt}) = \mathbf{0}
     * \f]
     *
     * As every DAE can be written as a first order DAE, it's acturally the description of a general
     * DAE system. Numerical approaches for the solution of DAEs can be divided into roughly two classes:
     * (i) direct discretizations of the given system and (ii) methods which involve a reformulation
     * (e.g. index reduction), combined with a discretization.
     *
     * The idea of a direct discretization is simple: approximate u and u′ by a discretization formula
     * like multistep methods or Runge-Kutta methods. The solution must be accomplished with a type
     * of iterative numerical method such as a Newton method. These technical difficulties are why,
     * in general, a direct discretization of fully implicit DAEs of index higher than one is not
     * recommended.
     *
     * It should also be noted that the implementation of direct discretization methods for DAEs
     * faces some additional practical difficulties such as how to obtain a consistent set of
     * initial conditions, a treatment of the ill-conditioning of iteration matrix, and finally,
     * error estimation and stepsize control for index-2 Hessenberg DAEs.
     *
     * @see http://www.scholarpedia.org/article/Differential-algebraic_equations
     *
     * PETSc’s default method for solving the nonlinear equation is Newton’s method. The general form
     * of the n-dimensional Newton’s method for solving \f$ \boldsymbol{F}\left ( \boldsymbol{x} \right )
     * =\boldsymbol{0} \f$
     *
     * \f[
     *  \boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-\boldsymbol{J}\left ( \boldsymbol{x}_{k}\right )^{-1}
     *  \boldsymbol{F}\left (\boldsymbol{x}_{k} \right )
     * \f]
     *
     * @see https://petsc.org/release/manual/snes/
     *
     * Thus, if we solve equations by Newton’s method, then we will need to differentiate F, either
     * exactly or by finite differences.  An implicit method will generally evaluate F a variable number
     * of times per time step, depending on the number of Newton steps
     *
     * @see Bueler E. PETSc for Partial Differential Equations: Numerical Solutions in C and Python[M]. 2020.
     *
     * In this form, however, only fully implicit solvers are appropriate. Try the IMEX (IMplicit-EXplicit)
     * methods to split the operators into two parts,
     *
     * \f[
     *  \mathbf{F} (t, u, v) = \mathbf{G} (t, u)
     *  v = \frac{du}{dt}
     * \f]
     *
     * where the functions F and G represent the implicit and explicit parts of the operator, respectively.
     * For “stiff” problems or those with multiple time scales F will be treated implicitly using a method
     * suitable for stiff problems and G will be treated explicitly when using an IMEX method like TSARKIMEX.
     *
     * @see https://petsc.org/release/manual/ts/
     * @see Ascher U M. Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. 
     * Siam Journal on Numerical Analysis, 1995, 32(3):797-823.
     *
     * @author nene
     * @date November, 2024.
     *
     */
    class ALGEBRA_EXPORT ODESystem
    {
    public:
        explicit ODESystem ();
        virtual ~ODESystem ();

        // // The number of equations in the system
        // virtual int GetEquationCount () const = 0;

        // /**
        //  * Take Template Pattern to maintain compatibilities with Boost.Odeint
        //  * @see https://www.boost.org/doc/libs/1_85_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts/system.html
        //  * @see GoF's Template Pattern
        //  */
        // void operator() (const Vec_t u, Vec_t dudt, Real_t t);

        /// The value of F( @a t, @a u, @a v) at the time @a t.
        virtual void ComputeImplicitValue (const Real_t& t, const Vec_t& u, const Vec_t& v, Vec_t& f) = 0;

        /// The value of G( @a t, @a u) at the time @a t.
        virtual void ComputeExplicitValue (const Real_t& t, const Vec_t& u, Vec_t& g) = 0;

        /// The jacobian of F( @a t, @a u, @a v) at the time @a t.
        virtual void ComputeImplicitJacobian (const Real_t& t, const Vec_t& u, const Vec_t& v, SpMat_t& dfdu, SpMat_t& dfdv) = 0;

        /// The jacobian of G( @a t, @a u) at the time @a t.
        virtual void ComputeExplicitJacobian (const Real_t& t, const Vec_t& u, SpMat_t& dgdu) = 0;
        /**
         * Solve the @a v that satisfies the algebraic equation F( @a t,  @a u, @a v) = G( @a t, @a u) at the time @a t.
         * For \f$ \dot{u} = G\left ( t,u \right ) \f$, it's the value of G( @a t, @a u) at the time @a t.
         */
        virtual void ComputeDerivatives (const Real_t& t, const Vec_t& u, Vec_t& v) = 0;
    };

    /**
     * @brief Represent the second-order ODEs \f$ \ddot{u} = G\left ( t, u, \dot{u} \right ) \f$.
     *
     * The general form of the the second-order ODEs can be described as followings,
     * \f[
     *  \mathbf{F} (t, u, v, a) = \mathbf{G} (t, u, v)
     *  v = \frac{du}{dt}
     *  a = \frac{d^2(u)}{dt^2}
     * \f]
     *
     * @author nene
     * @date November, 2024.
     */
    class ALGEBRA_EXPORT SecondOrderODESystem /*: public ODESystem*/
    {
    public:
        explicit SecondOrderODESystem ();
        virtual ~SecondOrderODESystem ();

        // // The number of equations in the system
        // virtual int GetEquationCount () const = 0;

        // /** @name Dedicated to the structural dynamics governing equations. */
        // ///@{
        // virtual void ComputeMassMatrix (SpMat_t& M) = 0;
        // virtual void ComputeLumpedMassMatrix (Vec_t& M) = 0;
        // virtual void ComputeDampingMatrix (SpMat_t& C)  = 0;
        // virtual void ComputeStiffnessMatrix (SpMat_t& K) = 0;
        // virtual void ComputeForce (Vec_t& F)             = 0;
        // ///@}

        // /**
        //  * Take Template Pattern to maintain compatibilities with Boost.Odeint
        //  * @see https://www.boost.org/doc/libs/1_85_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts/second_order_system.html
        //  * @see GoF's Template Pattern
        //  */
        // virtual void operator() (const Vec_t x, const Vec_t v, const Vec_t a, Real_t t);

        /// The value of F( @a t, @a u, @a v, @a a) at the time @a t.
        virtual void ComputeImplicitValue (const Real_t& t, const Vec_t& u, const Vec_t& v, const Vec_t& a, Vec_t& f) = 0;

        /// The value of G( @a t, @a u, @a v) at the time @a t.
        virtual void ComputeExplicitValue (const Real_t& t, const Vec_t& u, const Vec_t& v, Vec_t& g) = 0;

        /// The jacobian of F( @a t, @a u, @a v, @a a) at the time @a t.
        virtual void ComputeImplicitJacobian (const Real_t& t, const Vec_t& u, const Vec_t& v, const Vec_t& a, SpMat_t& dfdu, SpMat_t& dfdv, SpMat_t& dfda) = 0;

        /// The jacobian of G( @a t, @a u, @a v) at the time @a t.
        virtual void ComputeExplicitJacobian (const Real_t& t, const Vec_t& u, const Vec_t& v, SpMat_t& dgdu, SpMat_t& dgdv) = 0;
        /**
         * Solve the @a a that satisfies the algebraic equation F( @a t, @a u, @a v, @a a) = G( @a t, @a u, @a v) at the time @a t.
         * For \f$ \ddot{u} = G\left ( t,u,\dot{u} \right ) \f$, it's the value of G( @a t, @a u, @a v) at the time @a t.
         */
        virtual void ComputeDerivatives (const Real_t& t, const Vec_t& u, const Vec_t& v, Vec_t& a) = 0;
    };
    
    /** @} */
}  // namespace Algebra

SG_NAMESPACE_CLOSE

#endif  // Algebra_ODESystem_h